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ABSTRACT 


The one-speed, steady state, diffusion equation was solved for a point 
Squrceand for a normal pencil beam in a homogeneous, isotropically 
Scattering slab. Numerical results obtained using diffusion theory 
were compared to available transport theory results for two slab thick- 
nesses. 

This comparison demonstrated that the diffusion theory approximation 
to the transport equation will yield accurate results except within 
about one half mean free path of a boundary and except within about 
three mean free paths of a source. The best agreement between diffusion 
theory and transport theory is obtained if the first radial buckling 
constant in the diffusion solution is chosen equal to the radial buckling 


computed using transport theory. 
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Regions of Inaccuracy for Diffusion Theory 


IT. INTRODUCTION 


After being introduced into a diffusing medium, neutrons suffer 
numerous collisions with the atomic nuclei of the medium. The collision 
process is quite complicated, since either scattering, fission, or 
absorption can occur. As a result of streaming and repeated collisions, 
neutrons are constantly changing their location in the medium as well 
as their velocity (speed and direction). Indeed, some neutrons disappear 
(absorption) while others reappear (fission) while many simply change 
their speed and direction (scattering). The motion of any one neutron 
would appear to be quite random. 

The statistical concept of a distribution function is convenient 
for describing this complicated motion. In essence, we consider "typical" 
neutrons and try to find the neutron density throughout a medium by using 
the principles of transport theory. Although these principles are 
Straightforward and the exact equation (Boltzmann equation) governing 
transport phenomena can easily be derived, the solution to this equation 
is usually quite complicated. [1] 

Only in recent years have methods been developed to obtain exact 
solutions to the transport equation. These methods, while yielding exact 
results, are highly complex and require considerable computer effort. [2] 
For these reasons many problems in reactor physics and shielding design 
have been approached using the diffusion approximation to the transport 
equation. However, diffusion theory is expected to yield accurate results 
only when the assumptions of Fick's law apply to the problem. [3] 

Quite recently A. Leonard and G. Garrettson [4] developed techniques 


to solve the neutron transport equation exactly for the class of problems 


involving multidimensional neutron sources in a one-dimensional medium. 
Because some numerical data was available for a certain subclass of 
these problems, it was felt worth-while to develop a solution to the 
Same subclass of problems using the diffusion approximation and to 
compare the results with the transport solutions. 

In this thesis, the Green's function for neutron diffusion in a 
homogeneous, isotropically-scattering slab, surrounded by a non reflecting 
medium, was determined using the diffusion approximation to the transport 
equation. In Chapter II, an analytical expression was obtained for the 
neutron density from a point eounee arbitrarily located in the slab. This 
was integrated to yield the neutron density from a pencil beam normally 
incident to the slab. A computer program was written to evaluate these 
expressions, and some of its features are described in Chapter III. In 
Chapter IV, the numerical results from transport theory are compared 
with those of diffusion theory. 

This thesis has two objectives. First, it was hoped that by comparing 
the numerical results of diffusion theory with those of transport theory, 
accurate estimates could be made of the region in which the diffusion 
approximation can be expected to yield satisfactory results. Second, it 
was hoped a systematic method could be found to choose the parameters 
used in the diffusion equation to yield optimal agreement between the 


results of diffusion theory and transport theory. 


II. THE DIFFUSION THEORY SOLUTION 


A. THE EQUATION OF CONTINUITY 

Under the assumptions of one-speed, steady-state, and isotropic 
scattering, the equation of continuity in the diffusion approximation 
US 


Dv*o(r) - (Zz, - vEe)o(r) + S(r) = 0 2.1 


a 
where D is the diffusion coefficient; #(r) is the one speed neutron flux 
as a function of position, oe is the one-group absorption cross section, 
Ze is the one-group fission cross section, v is the average number of 
neutrons produced per fission, and S(r) represents an arbitrary source 
distribution function (those neutrons which have not yet suffered a 
collision). S(r) will later be considered a point source. Under the 
assumptions of Fick's law, the diffusion coefficient is given by 


D 


i} 


E/ 324, where E is the one group scattering cross section and 
SR ae is the total cross section [3]. 

The physical interpretation of eqn. (2.1) is as follows: 
-Dv*o(r)d?r represents the leakage rate out of a volume element, d°r, 


Via streaming; Eoo(r)d°r is the absorption rate in d°r3vz,e(r)d°r_ is 


‘i 
the neutron production rate in d°r from fission; and S(r)d°r is the 
production rate in d°r from any other source. 

To obtain the diffusion approximation to the transport equation one 


must make the following assumptions (Fick's Law): 


1) The medium is infinite 
) The medium is uniform 
) There are no neutron sources in the medium 


( 

(2 

(3 

(4) Scattering is isotropic in the laboratory coordinate system 
(5) The neutron flux is a slowly varying function of position . 
(6 


) The neutron flux is not a function of time. 
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The medium we consider is a homogeneous slab of thickness t < ~ 
with isotropic scattering, and we solve the problem under the assumptions 
of one speed and steady-state. Thus assumptions (2), (4), and (6) are 
Satisfied but (1), (3), and (5) are not. 

In order to obtain the diffusion equation (2.1), it was necessary to 
assume the diffusing medium to be infinite in all directions [3]. How- 
ever, neutrons arriving from more than a few m.f.p. (mean free paths) 
contribute little to the neutron current density at a given point, so 
the diffusion equation is expected to be reasonably valid except within 
a few m.f.p. of the boundaries. 

The no source assumption, implying that all neutrons contributing 
to the neutron density are the result of collisions, is certainly not 
valid in the cases studied. The slabs we considered contained highly 
singular sources: the point source and the normal pencil beam source. 
However, Since few neutrons originating from either source can be expect- 
ed to survive more than a few m.f.p. without suffering a randomizing 
collision, diffusion theory should produce satisfactory results at 
distances greater than a few m.f.p. from a source. 

The assumption of slowly varying flux is not satisfied near the 
sources or boundaries. In the derivation of Fick's Law, only first 
order terms in the Taylor expansion were retained [3]. Although second 
order terms will not contribute to the current density, terms containing 
third and higher order derivatives will make a contribution. Thus, it 
Ts necessary to restrict Fick's Law to regions where the second derivative 
of the flux does not change rapidly with distance. This specifically 
excludes regions near the sources we consider since both the normal 


pencil beam source and the point source (both represented using Dirac 


10 


delta functions) are highly singular. It should also be noted that 
since flux tends to vary rapidly in strongly absorbing media, it is 
also necessary to restrict Fick‘s Law to systems in which zr, << Z.. 
Such a restriction yields an average number of secondaries per collision, 
ee ND 
Coes - f 
close to one [1]. 
From the preceding discussion, we can see that the diffusion theory 
approximations should produce satisfactory results in regions of the 
slab not too near boundaries or sources. Hopefully, our numerical results 


will indicate more precisely the region in which diffusion theory is 


accurate and the degree of inaccuracy outside this region. 


B. THE BOUNDARY CONDITIONS FOR THE STEADY STATE DIFFUSION EQUATION 
As is normally done in reactor boundary problems, the boundary 


conditions at the surface were chosen such that 


do 


Lae. .t 2.2 
where oe is the normal derivative of the flux and 2 is the extrapolation 
length [3]. This is equivalent to saying that the flux vanishes at a 
distance 2 from the boundary. The value of % was chosen such that the 
solution to the diffusion equation (2.1) matched as closely as possible 
the more rigorous solution to the transport equation within the slab 
(away from the boundaries where the diffusion equation is valid). 

For a planar free surface, transport theory shows that 2 = 0.71 ev gives 
the best match, where der is the transport m.f.p. [3]. This choice of 


2 should provide gcod agreement between transport theory and diffusion 


theory in the interior of the slab. However, since the diffusion equation 


1 


is not valid near boundaries, a solution obtained by this device will 
not give the correct density near the boundaries (e.g., the flux does 
not really vanish at a distance % outside the surface). 

An additional boundary condition is obtained in a subcritical medium 
(c <1). The neutron flux decreases toward zero as the distance from 


the source increases. 


C. SOLUTION OF THE DIFFUSION EQUATION 

Consider the homogeneous slab of thickness, +t, which has a total 
cross section, £, and which emits c secondaries per collision. The slab, 
infinite in both transverse directions, iS surrounded by a vacuum or 
pure absorber. [See Fig. 1.] Under the assumptions of one-speed, steady 
state, and isotropic scattering, the equation of continuity in the dif- 
fusion approximation is given by (2.1). 

Since o(r), the one-speed neutron flux, is a function only of the 


neutron density, n(r), and the neutron velocity, v (assumed constant), 


o(r) = vn(r), 


ean. (2.1) can be written 


Dv? n(r) - (r. ~ vee )n(r) +7’ = 0. Lee 


Where xe(2,7 + 2), and y, z €(-~,~), 
Tt 1s convenient here to introduce cylindrical coordinates, 
(x,r); r= (r,¢), and r= Vy y2 + z2, which are illustrated in Fig. 2. 


Since it is assumed that the external source 


Sie) >+ 0, for Xe(2,t + g), oe(0,27) 


| fae 


the boundary conditions for our problem are 





PURE 
ABSORBER 


PURE 
ABSORBER 


FIGURE |. A BARE HOMOGENEOUS SLAB 


y 





FIGURE 2. CYLINDRICAL COORDINATES 
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n(O,r) = n(t + 22,r) = 0 for F c(0,~), de (0,2n) 2.4 


and 


n(x,P) > Q FOY Xmen +22 )5. 660,27) 2.5 


yo 
To obtain the neutron density n(x,r) from any uncollided source, 
S(x;r), it is convenient to have the Green's function, GUxa is ee) 
meveten here as a function of cylindrical coordinates [ref. Fig. 2]. 
Then the formal solution to (2.3) for any source S(r)/v is given by 
ee 
n(x,¥) = dx' dq' OMG aa 750 
0 
The Green's function has the following properties: 


(1) G satisfies the homogeneous differential equation 


2 -_ -_ — 
Dy2G (s. vE¢)G 0 Var | 
for all 


X e(L,7 + 2), y, z e(-~,~), and x # x', r # ae 


(2) G satisfies the boundary conditions 


G(2,r; x',r') = G(r + 2,43 x',r') = 0 2.8 


for all xe(z,t + 2) and y,z,y',z' e{-@,~); [5] and G(x,r' x',F')—> 0. 


Yoo 
(3) G satisfies the proper jump condition at r = r'. In cylindrical 
coordinates ,[9] 


G(x.r; x',r') > - log |F-F'] 2.9 
ae 2 ea 
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Equivalently, G satisfies the diffusion equation 
DAKE a ON vEg)G(r3r') otra) = 0. a2 ..10 


with a point source, S(r) = 63(r - r'), where 


Mum Ges) 2: aere ete) cand "4.2 6 (-cie) - 
63 (r-r') is the three-dimensional Dirac delta function given by 


63(r - r') = 6(x - X') 6(y -y') 6(z - z') in cartesian coordinates and 
63(r - r') = 62(r - r') 6(x - x') in cylindrical coordinates, where 


~~ 
t 


gy elle SR ee 


a ae 
sae = 


[2 


To solve for the Green's function one can find the eigenfunctions to 
the homogeneous form of equation (2.1) by separation of variables. Then 
G(r; r') can be expanded in terms of these eigenfunctions, and the 
expansion coefficients can be chosen so that the jump condition is 
satisfied [7]. 


In .cylindrical coordinates equation (2.10) is 


Assuming a homogenecus medium and a point source located at x' on the x 
axis, symmetry dictates that G(x,r; x',0) is a function of r but not of 
¢. In general then, G is a function of |r - r'| rather than ¥ and T' 


separately. Defining 


g(x, | 


[st 


Stee heen a eGlbeis x 


we obtain 
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2 ~ re 2 
+ =| ieee ck) Chg efor -74 (r) S(x -x') 2.12 





sj — 
Q? 
Ne 
=) 
Slee 


where we have defined the coefficient 


> = vet 
Q BoE 2.13 


itl 


g(x..r3x') satisfies the boundary conditions (2.8). For fixed x‘e(0,7) 
and 7'(0,~), the homogeneous form of equation (2.11), satisfied by the 


eigenfunction w(x.r), is 


zy 5 2 mt a 
{2a + oe + xr | w(x ,r) = Q w(x.r). 2.14 


To find the eigenfunction p(x,r) we separate variables 


v(x sr) = X(x)R(r), 
and substitute into equation (2.14). This leads to the following 


separated eigenfunction problems: 


d2X L2X = 
ane : ZA 
with boundary conditions 
X(0) = X(t + 22) = 
and 
dR ,1dRy op - 
az x de + a“R = 0 2.16 


with the boundary condition 





* Note: For simplicity, the solution is forced to o at x = 0 and 
X ='1 + 22, rather than at t.= - 2 and x = tte. For numerical work 
and comparison with transport theory, a change of variables, Se aR. 


is introduced where X ¢(0,1t) iS inside the slab. 
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where 
a2 = Q + Lé. 


Equation (2.15) has solutions of the form 


Xx ) = C sin(L_ x) + C cos (L_x) 


where the boundary conditions require that C,, = 0 and that the eigen- 
values are 


eT: a aoteae 
Ls Reap: Mm =sle2,3; 


: _ 5 
By changing variables such that p=ayr, where Orn gf (he Ls 


equation (2.16) can be reduced to the modified Bessel's equation 


p2 d2R(p) dR(p) 


Se tp Gye - p? Rip) = 0. 2 Vi 


Since the solutions to equation (2.17) are the zeroth order modified 


Bessel functions 1 (p) and K (Pp); [6] the general solution is 


R(p) = C41. (p) + CyK, (p). 
The boundary condition requires that C, = 0 since PO 00 
De 


Expanding the Green's function in terms of the eigenfunctions 


HI 





XX) Sul ( o. 


and R_(¥) = Kar) 


we obtain 


ar A(x') X. (x) RACBE = feo lhe 2.19 
Since 
X (x) Ror) = v_(xs¥) 
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satisfies the homogeneous equation (2.14) and the boundary conditions 
(2.8), the expansion (2.19) for abies es) aly satisfies the 
conditions (2.7) and (2.8) required of the Green's function. It 

remains only to be shown that (2.19) satisfies the jump condition (2.3). 
This is obtained immediately since [6] 

ie, |ee |) = en |r - r’| 


So that R tr) satisfies the following differential equation [8]: 


aR, 
Sie 


s(— 


dR (1) 
spo 
Thus equation (2.19) is the proper form for the solution. 


Substitution (2.19) into equation (2.10), and using the Fourier sine 


representation [8], 
ad 


S(x -x') = a. (me, | sin( ,, 2.21 


m=} 


as well as equation (2.20), yields 


Py 


1 = ] . be 
A(x ) = Dit + 22 sin (DE), 


where we have utilized the orthogonality property 





Table 
: 0 form7n 


; max . Nn1rXx 
sin ( ! sin | dx = 
t + 22 t + 22 = 2b OR 


0 








form=on. 


Therefore, recalling that G(x,r; x',r') = g(x,|r - ¥']5x'), 
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MmX 


aD 
g(x, |r-r']sx") = = sin (SE } sin (SR) Ky (ag BEN) 
qe | . 
2.22 


is the required Green's function for the point source problem. 

It would, of course, be possible to develop the Green's function 
for the normal beam problems in a like manner. However, this is not 
necessary since equations (2.22) and (2.6) represent the formal solution 


to (2.3) for any source. In particular, 
aie 


o(x,%3x!) = g(x,|F-F']; x") 62(F)e * ax! Farag! 
O07 £ 2403 


is the neutron density in a slab from a pencil beam normally incident to 


the slab at x = y = z = 0, and this function satisfies equation (2.3) 


< Sige) ~ = ' me a 
with aa 6*(rje BS Tins o(x,|r-r'],x') is the Green's function 
for the class of problems involving neutron beams normally incident to 
a slab. Substitution of (2.22) into (2.23) and subsequent integration 


yields 





00 
o(x,r3x') = 7 Peete Say sin ( Tz} K (oF). 2.24 
m= 

To facilitate comparison of our results with those of transport 
theory [1], we measure length in dimensionless units of mean free path 
(1 m.f.p. = 1/z) and shift coordinates in x and x', x > x - 2, So that 
x e(0,t) lies inside the slab. Then the solutions (2.22) and (2.24) 


take the form 


CO 
\ 
wt) on mmx 8) ) .  ( ma(xt2) y 
qg(x,r3x') =C / sin ( om te) sin ( rig y ) K (oF) 2.25 
Ma! 
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and 


OO 
) -(t+22) | 
we = \ ] = (-] )Me (r+ . m x+2) ~ 
BX TSX ) = (C tat eens Sun [ mele) | K (ot) » Zao 
m= 


where all lengths are expressed in units of mean free path. For example, 


2 
a mir ~ 
a = \/ tr + (a) (rr). 2.27 


where rr = r/(1/z) is the radial distance expressed in mean free paths. 
For equations (2.26)-(2.27), and AAG this point on, all distances are 
understood to be expressed in units of mean free path (e.g., £(1t+22) > 
t+22 m.f.p.). The constants C and C' will be adjusted to normalize the 


diffusion theory results to transport theory results at some chosen point. 
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IIL. NUMERICAL TECHNIQUES 


A Fortran IV program (appendix 2) was written to compute numerical 
results using equations (2.25) and (2.26) for comparison with numerical 
data obtained by transport theory. Since both equations involve an 
infinite series, it was necessary to truncate the series at some point 
riborder to obtain a solution. The purpose of this chapter is to 
justify the criteria used to terminate the series and also to explain 


the choice of the particular values used for the parameter By orn 


A. TRUNCATION OF THE SERIES 

Both equation (2.25) and (2.26) are infinite series solutions. In 
order to obtain numerical output, it was necessary to truncate these 
series at some point in the summation process. That is, 


ae) N-! 


i x Te 
Lo Aan PK (on MY (Aan OK olen) = Sy sal 
where we choose N large enough so that Ry/Syy is less than some chosen 
e, where 
Ry = A (x)K, (ant) 3.2 


is the remainder term. ine maximum absolute value of the A_X (x) terms 
mom 


in both equation (2.25) and (2.26) is unity. Thus the Ky Bessel function 


term dominates the series for we >> 1. 


since 
00 
AnXml 1 Ry s YK lant) 

ake Ue ns ys, oh n® eae: 

| 7 M=N 

factoring out K anh) and noting that [6] 
‘ Fae 
K(p) ay 3.4 
Doo 


2] 


with 


[itr c 


yields 





7 “es k 
my, Tesi (Tan) | 
Ry <= K (ayn) KOEN C / ; 3.4 
K= 0 


This power series can be summed to yield a convenient estimate for an 
upper bound of the remainder: 
ae K loyr) - | 
Nes ee le on0 
1 -e ‘+22 
In the computer program the series is truncated at N terms with N 
large enough so that 
—— cant S210 
where we chose ec sufficiently small to give the desired accuracy. Then 
(3.5) yields the upper bound 
Ry = 


SS 


Sie Sey 
SN 1 - exp[-2,,] 


for the remaining fraction. This upper bound is of the order of e« except 





when r << wes , so that the choice c« = 0.0005 should give at least 

three place accuracy except for small Yr. However, (3.7) represents a 

very conservative estimate of the upper bound, so we would expect good 
accuracy with this choice of ce even when r is small. In fact; for Y = 0.05 


and N large enough so that 


K (on) 
N 


< Q.0005, 


22 


R 
numerical results indicate that es < 0.001, while (3.7) yields the 
N 


R: 
Conservative upper bound < (Ole 
N 


Therefore, on the basis of this discussion, it would appear that 
(3.6), with « = 0.0005, is a sufficient truncation requirement to insure 


three-place accuracy. 


B. CHOICE OF THE PARAMETER Q/z4 

In equations (2.25) - (2.27) the parameter Q/z2 appears, where Q is 
defined in eqn. (2.13). This parameter is a function of the one group 
material cross section of the slab. However, the one-group cross sections 
depend specifically upon the energy-dependent flux, which is not avail- 
able. Therefore, one would like to choose Q/z2 so that diffusion theory 
will yield results that agree as closely as possible with the results of 
Oone-speed transport theory. In this thesis it is demonstrated how Q/:2 
might be chosen to accomplish this. 

For comparison, this parameter was chosen in two ways. The first 
choice assumes diffusion theory for a medium in which there is no fission. 
For the second choice, it is noted that as * increases, the transport 
theory solutions and the diffusion theory solutions decay like K (Bt) 
and K (ar); respectively, where 8, are the transport theory radial 
buckling modes [2]. To force diffusion theory to behave like transport 
theory for large r, the first buckling modes in diffusion theory and 
transport theory are chosen equal, Oy = 81° 

To obtain the expression for Q/s4 in the first case, recall the 
parameters 


eer >: 
g= 2 ft 3.8 


and 


(ee: 


Se i a 3.9 


defined in Chapter II. In the diffusion approximation [3] D = Eo /3r? : 


and if there is no fission (Ze =D). Q/z2 reduces to 


: a _ 3(1-c) 3.10 
C : 


since c= £./(Z, + £,) in this case. 

Recall that the diffusion approximation demands a highly scattering 
medium, c 1.0, so c = 0.9 is chosen for the numerical work presented 
here. For this value of c, eqn. (3.10) yields the value Q/z2 = 0.333. 

For the second case the parameter Q/£* is obtained by comparing the 
transport theory and diffusion theory solutions. For large transverse 
distances from the source, lr-r' | >> 1, a relatively simple transport 
theory solution is obtained [2] since the finite (M < ~) point spectrum 


modes dominate the continuous spectrum modes: 


ba 
K (g_lr-r' ee 
G, .(rsr') S Pionuiies r(x) r(x) oC e FEET} ain 


/ 


i 

where r(x); m= 1,2,+:+,M, are the transport eigenfunctions for the 
discrete eigenvalues 0 < 6) < By < *'* < By <1 [2]. Using eqn. (3.4) 
it can be seen that the first mode dominates (3.11) for large radial 


arguments, |r-r'| >> 1: 
-85|r-r'| 
re Cuureesa ens i e 
mets ) si K (8, [r-r' |) ry (x) ry (x ) 80 | Sel 
fifee, 
Likewise, the first mode m = 1 dominates the diffusion theory solution 
(2.25) for |r-r'| >>1. Therefore, if we choose the first radial buckling 
modes to be equal, a, = By> the two solutions will behave similarly for 


NS ff 
ieee role>> 1. 
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ieee. Q ; tal . 
With this choice for 4 = } 3 a5 {| , the following 


expression is obtained for Q/2?: 


_ 2 
oy = Bs (<5 } | 3.12 


Numerically this corresponds to Q/r4 = 0.2757, for c = 0.9 and t = 10 
and to Q/z2 = 0.1918, for c = 0.9 and t = 3. [2] 
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TV. DISCUSSION OF RESULTS AND CONCLUSIONS 


The following results were obtained using the computer program discus- 
sed in Chapter III. The cases presented are those for which the exact 
solution has been obtained by G. Garrettson [2] using transport theory. 
The purpose of this chapter is to compare the results obtained using dif- 
Bein theory with those obtained using transport theory and to evaluate 
the performance of diffusion theory for problems of this class. These 
comparisons were done for saeRe of homogeneous slabs of thickness 1 = 10 
and t = 3m.f.p. (mean free paths) surrounded by a non-reflecting medium, 
under the assumption of one speed, steady state, and isotropic scatter- 
ing. All graphical results (Figs. 3 to 17) discussed in this chapter 
are found in Appendix 1. 

The first results presented, Figs. 3-6, are those for a pencil beam, 
normally incident at (0,0,0) to a slab of thickness + = 10 m.f.p., with 
a multiplication constant c = 0.9. The data has been normalized to the 
transport theory results at r = 4.4448, x = 5.1282 in order to facilitate 
the comparison. Figures 7-10 depict the results obtained from a pencil 
beam normally incident at (0,0,0) to a slab of thickness t = 3 m.f.p. 
with c = 0.9. Figures 7 and 8 were normalized at r = 2.01, x = 1.5, 
while Figs. 9 and 10 were normalized at r = 0.2, x = 1.5. Finally, Figs. 
11-16 represent a comparison of the point source solutions in a slab of 
thickness t = 10, with c = 0.9, for point sources located at depths of 
1 and then 5 m.f.p. 

In all but Figs. 11-13 the diffusion solutions are presented for both 
values of Q/z4 discussed in Chapter III. In this chapter, and in Appendix 


1, Q/z4 is referred to simply as Q, since 5 = 1 in units of mean free path. 
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A. COMPARISON OF NORMAL BEAM SOLUTIONS 

It is interesting to note the evolution of the neutron density from 
a pencil beam, N(x,r), vs. x e[0,1], as the radial distance, YT, increases. 
For very small r (ref. Figs. 3, 9, and 10), the neutron density tends to 
have the same shape as the once-collided neutron density, which in turn 
resembles the uncollided density (a step function times an exponential). 
Transport theory predicts a very sharp discontinuity near the boundaries 
which is physically explained by the surface leakage [2]. Although dif- 
fusion theory correctly yields an exponentially decaying density within 
the slab, it does not yield the steep gradient given by transport tneory 
within about one half m.f.p. of either surface. The results of the two 
Brennics diverge partly because transport theory accounts for the 
preferential streaming toward the boundaries while diffusion theory does 
not. In fact, as previously discussed, the diffusion equation is simply 
not valid near boundaries. (Ref. Chapter I1). 

For any x e(0,1), the diffusion theory solution appears to be most 
inaccurate for small radial distances, ft <1. This inaccuracy is apparent 
from both the normal beam and point source data (ref. Figs. 3, 9, 10, 11, 
and 14), and it arises partly because diffusion theory does not consider 
the preferential streaming from a singular source in the region of the 
source. Since Fick's Law is derived under the assumption that there are 
no sources in the medium (ref. Chapter II) we should not expect the 
diffusion results to be accurate near sources. In fact, within one m.f.p. 
of the source, the errors induced by using diffusion theory may be as 
large as 80 percent. However, beyond about three m.f.p. the two theories 


yield nearly identical results (if Q is computed using eqn. (3.12)). 


* Recall that a straight line on a semilog plot corresponds to an 


exponential. 
2/ 


B. COMPARISON OF POINT SOURCE SOLUTIONS 

For solutions to the point source problem, the most noticeable 
discrepancy between the transport and diffusion data appears in the 
region within about one m.f.p. of the source (ref. Figs. 11 and 14). 
As previously mentioned, diffusion theory's inaccuracy near sources 
is partially due to its failure to consider the preferential streaming 
in -the region of the source. As the distance from the source is 
increased, the relative error decreases unless a boundary region is 
encountered. As expected, within one half m.f.p. of the boundary, the 
diffusion theory solution begins to diverge until a discrepancy of 


about 20 percent is noted at the edges. 


C. CHOICE OF THE PARAMETER Q 
It should be noted (ref. Figs. 5, 7, 8, and 16) that for large 
radial distances, r, the agreement between transport and diffusion 
theory can be greatly enhanced by the proper choice of the parameter 
3(1-c) 


Q (ref. Chapter III). The diffusion theory yields Q = oral for a 


medium in which there is no fission. However, it was expected that 


2 
= g2 . (_ 7 
Q= ot - (5 iz 
where By is the first transport eigenvalue [2], would yield better 
results. By examining both the normal beam and the point source data 
for the two choices of Q, the most optimum choice is apparently given 
by equation (3.12). This choice forces the diffusion solution to decay 
in Y at the same rate as the transport solution (for large r) because 
the first modes dominate for >> 1. The numerical data bears this out 
Tt 
TH2R 
solution for t > 3. Unfortunately, this choice of Q does not consistently 


Since the Q = By -[a3° solution is nearly identical to the transport 
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improve the results near the source (or near the boundaries) since the 
assumptions for Fick's law are not valid in these regions and since the 
first mode does not dominate the solution for small r. A plot of 8) VS 1 
is given in G. Garrettson's thesis for c = 0.9, and data for other c's 

is contained in an appendix [2]. 

A comparison of results in the t = 10 and 1 = 3 slabs (refs. Figs. 
3-10), indicates that diffusion theory is more accurate in thick slabs. 
This is to be expected since Fick's law is not valid near the boundaries, 
and in a thin slab the neutron is never very far from one of the two 


boundaries. 


D. CONCLUSIONS 

In general, one can say that the neutron density computed using 
diffusion theory can be expected to yield reasonably accurate results 
except as indicated in the shaded areas of Fig. 17. If one is not 
interested in the density within about three m.f.p. from the source and 
within about one-half m.f.p. from the edges, diffusion theory can success- 
fully be used. However, if one is interested in the density within these 
regions (e.g., the shaded areas in Fig. 17), it will be necessary to 
resort to the more rigorous and time-consuming techniques of transport 
theory to obtain accurate results. 

One shouod also consider that the diffusion theory solution reouires 
very long to compute for small radial distances from the source, since 
the smaller |r-r'| is, the more terms the truncated series includes. For 
example, computer runs which included r = 0.01 required as long as 18.5 
minutes to complete, which is comparable to the computation time required 


for the most difficult transport theory solutions. These runs, even 
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though lengthy, did not provide accurate data, and one must conclude 
that only the transport theory solution should be used to compute 
densities at very small radial distances. 

Lastly, the choice of the parameter Q is an important factor in 
obtaining agreement between the results of diffusion theory and transport 
theory in the regions where the diffusion equation is valid. The choice 
(3.12) significantly decreases the discrepancy between the two results: 
at large radial distances from the source, e.g., r > 3 m.f.p. 

In shielding problems with beam sources, the effect of these 
discrepancies will probably be most noted when calculating the transmis- 
Sion or reflection from the region of the beam. In these regions the 
diffusion results are most inaccurate since the emerging current from 
a region of the slab is a function of the neutron density in that region 
[3]. Diffusion theory should give satisfactory results for both trans- 


mission and reflection from regions not too near the beam. 


E. SUGGESTIONS FOR POSSIBLE EXTENSIONS 

For slabs sufficiently thin, no discrete poles exist for the transport 
equation [2]. In these cases, to choose a proper value for Q, an effec- 
tive mode of the transport equation would have to be found empirically. 
This could be done by making a semilog plot of the transport solution 
for fixed X, aS a function of r (for large r), and using the slope to 
determine the effective decay constant. Mathematically, this would 
correspond to the smallest "effective" pole in the continuous spectrum 
of the transport operator. These effective poles should be investigated, 


not only empirically, as mentioned above, but also mathematically. 
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Further, it would be interesting to investigate the results obtained 
over a wide range of values of both cand c. One would expect diffusion 


theory to be less accurate for small +t, because of leakage, and also for 


tok ee 
small c = oe which corresponds to a highly absorbing slab. The 


diffusion approximation is not expected to be valid unless ZR, << 2. 
(Ref. Chapter II). 
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APPENDIX 1 - Numerical Results 


Appendix 1 contains graphical displays comparing numerical data 
obtained using transport theory with that obtained using diffusion 
theory. Figures 3-10 represent semi-log plots of the neutron density 
in a slab, as a function of normal depth, from a pencil beam normally 
incident to the slab, for various radial distances from the beam path. 
Figures 11-16 represent semi-log plots of the neutron density ina 
Slab, as a function of normal.depth, from a point source within the 
Slab, for various radial distances from the source. 


This data is discussed in Chapter IV. 
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APPENDIX 2 - Fortran IV Program Listing 
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port theory results for two slab thicknesses. : 


This comparison demonstrated that the diffusion theory approximation to the 
transport equation will yield accurate results except within about one-half mean 
free path of a boundary and except within about three mean free paths of a source. 
The best agreement between diffusion theory and transport theory is obtained if 
the first radial buckling constant in the diffusion solution is chosen equal to 
the radial buckling computed using transport theory. 
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